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I. INTRODUCTION 

Quantum computation is the investigation of the properties of devices which use quantum mechanics to process 
information [1, 2]. Such quantum computers exist in the laboratory now, but currently are only able to process a few 
bits of information. Many experimental proposals for quantum computing exist - all of which share the goal of the 
controlled manipulation of the quantum state of the computer. These include superconducting systems [3], trapped 
atoms and ions [4, 5], Nuclear Magnetic Resonance (NMR) [6], optical implementations [7, 8], Rydberg atoms [9], 
quantum dots [10, 11], electrons on Helium [12] and other solid state approches [13, 14]. For certain specialized 
applications it is believed that quantum computers are more powerful than their classical counterparts [15-19]. 

One of the most promising areas of application for future quantum computers is quantum simulation: the use of a 
quantum computer to simulate or emulate another quantum system. Feynman's idea that quantum machines could 
be useful for simulating quantum systems was one of the founding ideas of quantum computing [20] . In the case of 
lattice Hamiltonians, specific techniques were first developed by Lloyd [21]. Zalka and Wiesner [22, 23] developed 
methods for wave mechanics and Lidar and Wang used this to define an algorithm for the calculation of the thermal 
rate constant [24]. Work to date in the area of quantum simulation has considered both static [17] and dynamic 
properties of quantum systems [22, 23]. Recently these ideas have been applied to the quantum statics and dynamics 
of molecular systems [25, 26]. 

Accurate computation of the ground state energy of molecules is a basic goal in chemistry, but an area where 
accurate results remain expensive and elusive in many areas. This is also an area where quantum computers are 
known to offer significant advantages over classical machines for which the cost of numerically exact calculation grows 
exponentially with the size of the problem [25]. In the context of physical chemistry quantum computers are suited 
to full-CI electronic structure calculations and to the simulation of chemical reactions without the Born-Oppenheimer 
approximation [25, 26]. In applicable cases, these quantum algorithms provide an exponential performance advantage 
over classical methods. Work to date indicates that a quantum computer with a few hundred qubits would be a 
revolutionary tool for quantum chemistry. This number of qubits is a realistic target for sustained experimental 
development. However, given the small number of qubits currently available any application of quantum computing 
faces the challenge - how do we get there from here? 

Quantum computation is slowly moving out of its infancy and it is now possible to consider the implementation of 
the smallest realizations of many of the basic algorithms known to present an advantage over classical computation. 
Significant experimental progress has been made in NMR [27] , in trapped atoms and ions [28-32] , superconductors [33- 
35], optical implementations of the gate model [36-39] and other approaches [40]. While the field may be maturing 
beyond its earliest days, there remains a large amount of work to be done before quantum simulation can compete 
with existing classical computers. How should we think about the early efforts to implement quantum algorithms for 
chemistry, and how much progress do they represent towards the goal of real quantum chemical discovery on quantum 
devices? This chapter proposes a framework for evaluating partial progress towards useful quantum simulators for 
quantum chemistry. The basic idea is to use the over sixty years of development of quantum chemical methods, 
implementations and results on classical machines to guide and evaluate the development of methods for quantum 
computers. 

A great deal of effort was focussed on quantum chemistry over the last sixty years, and it is possible to go back to 
the literature and find the calculations performed on early computers - some from the vacuum tube era [41]. These 
calculations are trivial by modern standards, but provide exemplars which could be implemented on early quantum 
computers (with, say, five to twenty quantum bits). By starting at this point, one can then trace the evolution of 
the classical calculations to the present day, and thereby lay out a roadmap leading from calculations which can be 
performed on quantum computers now to calculations at the research frontier. 

For calculations in which 1 — 20 qubits are used to represent the wavefunction we propose a historical approach 
to the development of experimental realizations of quantum chemistry on quantum computers. The first algorithm 
in this sequence has been implemented experimentally in linear optics quantum computing for a minimal basis set 
representation of the Hydrogen molecule [42, 43]. 

We first address the merit of such an approach: Why not simply plug in appropriate minimal bases from modern 
electronic structure methods for simple systems? One motivation is to reexamine the problems of quantum chemistry 
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from a new viewpoint. In spite of the polynomial scaling of these algorithms with the problem size on quantum 
computers, gate and qubit estimates for simulation remain large. Theoretical improvements to these algorithms could 
significantly reduce the experimental requirements for simulation of technologically interesting molecules. Methods 
adapted for classical computers may not be the best for quantum computers, and so looking back at methods which 
were perhaps abandoned is worthwhile. Another motivation is to avoid the temptation to focus on those problems 
that are particularly suited to quantum simulation, and to instead evaluate the methods against a pre-selected set 
of benchmarks. As such, these historical benchmarks provide a complement to valuable investigations of a set of 
standard chemical examples, as is performed in [44]. Finally, it is interesting to measure progress in the quantum 
implementation against the historical development of the classical methods. This gives a straightforward way of 
estimating the distance yet to be covered before quantum techniques can compete with classical methods. 

This example-by-example approach should be contrasted with more general theoretical inquiries. Much theoretical 
work in quantum information is focussed on quantum complexity theory [45, 46]. In particular, the Local Hamiltonian 
problem has been extensively studied [1, 47, 48]. The Local Hamiltonian Problem formulates the task of finding the 
ground state of a Hamiltonian in the form of a decision problem. As such, the problem can be studied in terms of 
two solution concepts. First, one may consider versions of the problem where the solution may be found efficiently. 
Classically such examples would fall into the complexity classes P (Polynomial time) or, if a probabilistic algorithm is 
required, they would fall into BPP (Bounded-error Probabilistic Polynomial time). Secondly, one may consider those 
problems for which a given solution may be verified efficiently. Classically such problems would fall into the class NP 
(Non-Deterministic Polynomial time), or again, if the verification procedure is polynomial, into MA (Merlin- Arthur 
- see below for an explanation of this nomenclature) . Of course, the question of whether the two solution concepts - 
finding versus verifying a solution - are in fact the same is the question of whether P=NP [49]. 

In the classical world it is known that there are Hamiltonians for which verification of the ground state energy 
is as hard as any problem in NP [50]. At least one classical algorithm for fermionic quantum systems, Quantum 
Monte Carlo, is also as difficult as any problem in NP in the worst case [51]. On the quantum side, it is known 
that the eigenvalue decision problem in general for a local Hamiltonian on qubits is as hard as any problem in 
QMA [1, 47, 48]. However, these complexity-theoretic results are mostly for worst-case complexity, and hence one 
may state the advantage of quantum algorithms for ground state properties in complexity theoretic terms as dependent 
on the existence of a subset of physically interesting problems that are quantumly tractable but classically hard. 

In the present paper we do not propose an approach to obtaining complexity-theoretical results on the efficiency 
of quantum algorithms for ground state problems. Such a project might be described as Newtonian or Cartesian - 
the determination of grand theory covering all possible cases. Needless to say, such lofty goals arc; difficult to acheive 
in practice. Complexity theoretic proofs of the advantage of many widely used classical algorithms are few and far 
between. Indeed, in the case of Density Functional Theory and Quantum Monte Carlo the proofs are in the other 
direction, showing that the worst cases of these algorithms are likely to be hard, just as is the case for the local 
Hamiltonian problem [1, 47, 48, 51, 52]. 

We may compare the Cartesian approach of quantum complexity theory to the approach of the present paper. We 
propose instead a Baconian approach, where specific examples are investigated in detail one at a time. These two 
approaches were recently contrasted by Dyson [53] , who characterized them in terms of Cartesian birds and Baconian 
frogs. The Cartesian birds fly over the landscape describing its broad features and general properties, whereas the 
Baconian frogs are happy to sit in the mud and investigate the details of what lies near them. Of course, we hope 
that these approaches are complementary What we learn from the detailed investigation of a set of specific examples 
should be of use in correctly defining the right conditions under which quantum algorithms for ground state problems 
can be rigorously proved to be efficient. Setting theoretical considerations to one side we also hope that the elucidation 
of a well defined set of examples, and their use as a yardstick of progress, will help in the further development of 
experimental procedures and hardware for quantum simulation. 

II. QUANTUM COMPUTING 

In this section we briefly review some salient techniques of quantum computation. Here we shall focus on brute 
force methods aimed at producing the minimal number of qubits - suitable for the first experimental implementations. 
Such implementations are useful for several reasons, not least because as they determine whether the experiments 
have sufficient precision to obtain the eigenvalues to a reasonable accuracy. For further details an the fully scalable 
approach we refer the reader to [25] and to the Chapter in this volume by Aspuru-Guzik et al, and [54, 55]. In those 
works the fully scalable approach to quantum simulation is reviewed in detail. 
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FIG. 1: A schematic circuit for one step of the phase estimation algorithm. The upper line is the readout register (taken to be 
one qubit here) and the lower line is the state register which carries an encoding of an eigenstate of the time evolution operator 
U. 



Phase estimation 



Phase estimation is a fundamental primitive of quantum computation [1, 2, 19]. Here we review those elements of 
the procedure that are relevant for our subsequent discussion. The circuit we shall analyze is shown in Figure 1. It 
is composed of two registers - the readout register - which we take to be one qubit, and the state register, which we 
consider to be an arbitrary number of qubits. The input state to the readout register is simply logical zero, |0), while 
the input state on the state register is an eigenstate of the Hamiltonian, and therefore of the time evolution operator: 



U\^)=e~ i * 



where 
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for evolution time r and eigenenergy E. A Hadamard gate (which is the Fourier transform over the group Z/2Z) is 
applied to the readout register, resulting in state: 
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The next stage is a controlled application of the unitary operator U corresponding to time evolution for a period r, 
resulting in state: 
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The final step before measurement is to again apply the Hadamard gate to the readout qubit, resulting in the state 

1 
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(5) 



It is worth noting that in the case where \ip) is exactly equal to an eigenvector of U, the readout and state registers 
are not entangled at all in states |ST), \S2), \S3). A measurement performed on the readout qubit when the system 
is in state 15*3) results in zero or one with the following probabilities: 



Prob(O) 
Prob(l) 
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1 + cos(^/2) 
1 - cos(0/2) 



(6) 
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We can now proceed from our description of phase estimation to energy measurement, following Abrams and 
Lloyd [56]. Suppose that we can estimate arccos((cr z )/2) = <j> to fixed precision. Let U be the time evolution operator 
of a system for time t: 

U = exp(-itH/K) (j) = tE/h 
suppose the energy scale is such that E/Ti < 1 so that we can write E/Ti as a binary fraction: 

E/n = 0.E 2 E iE8 E 16 ... = ^ + f + f + ^ + ... 

Z 4 o lu 

then choose i = 2n2 n so that: 

= iV27r + 7r£; 2n /2 + ... (7) 

Using the recursive phase estimation procedure define in [25, 54, 55] we can repeat this calculation to estimate E one 
bit at a time. This is the procedure that has been implemented for calculations of the Hydrogen molecule in [42, 43]. 



B. Time evolution and the Cartan Decomposition 

Evidently, the procedure outlined above requires the implementation of controlled operators corresponding to time 
evolution for a sequence of times 2V for i = 1,2,.... What is the difficulty of this problem if we know nothing 
about the form of U except its dimension, i.e. all that is known is the number of qubits in the state register? The 
problem we wish to address is the decomposition of an arbitrary unitary operator into a sequence of one and two 
qubit gates. These one and two qubit gates are the basic primitive operations available to the experimentalist, and so 
we wish to proceed from an arbitrary unitary to a quantum circuit over a particular gate set. A typical choice is the 
controlled not (CNOT) operation combined with arbitrary operations on single qubits - and this choice is known to 
be universal, i.e. sufficient to implement any unitary operator [57]. Much work in quantum computing has addressed 
this problem, including approaches that set bounds on the number of gates required [58, 59], and approaches which 
can constructively obtain the full circuit, including all parameters for the one and two qubit gates [57, 60-62]. For 
factoring an q qubit unitary over CNOT gates and arbitrary one qubit rotations the best known performance is: 

23 4" - 3 2« + 4 (8) 
48 2 3 W 

CNOT gates and 4 q — 1 one qubit gates, where a one qubit gate is either an X, Y or Z rotation on a single qubit. 

The number of gates for an arbitrary unitary grows exponentially with the number of qubits. This is expected for 
a method that pays no attention to any special structure that may be present in U beyond the number of qubits on 
which it acts and the tensor product structure of the Hilbert space of a set of qubits. Scalable implementations of time 
evolution operators must exploit the structure of the Hamiltonian in order to realize a polynomially scaling number of 
gates. Fortunately, for quantum chemical Hamiltonians in the second quantized formalism it is well understood how 
to do this, with a number of gates rising as q 5 . For details of the scalable approach see the article by Aspuru-Guzik 
in this volume, or [54, 55]. 

The bounds given in 8 arises from a technique referred to as the quantum Shannon decomposition [59]. The 
underlying mathematical tool used in much of the literature is the Cartan Decomposition of a Lie Algebra [63, 64], 
and the corresponding decomposition for the Lie Group - in this case the unitary group. 

What is the Cartan decomposition for a Lie group? Let us restrict to unitary operators of dimension U(2 n ) with 
determinant one, meaning the Lie group SU(2 q ). Let us call this group G, with the corresponding algebra being 
the antihermitian matrices (which physically correspond to Hermitian observables, times i). The first step in the 
decomposition is to identify a subalgebra, and hence a subgroup of G. We call this subgroup K. For our purposes we 
may consider two types of subalgebra (and hence subgroup). Firstly we can choose the subgroup to be the real unitary 
matrices, in other words the orthogonal group SO(2 q ) 1 with the corresponding subalgebra being real antisymmetric 
matrices. Secondly we can consider block diagonal unitaries with blocks of size p and r such that p + r = 2 n . It is 
natural for systems composed of qubits to choose p = r = 2 q ~ 1 , so that the blocks are half the size of the original 
blocks. 

Having selected our subalgebra we have the algebra g corresponding to the original group G = SU(2 n ), the 
subalgebra k with corresponding subgroup K. The orthogonal complement of k in g is called m, and for a Cartan 
decomposition we find that the commutators obey the following relations. Firstly, A: is a subalgebra, so that: 



[fci,A;2] £ k for k\ and k 2 in k. 



(9) 
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FIG. 2: The two steps in the Quantum Shannon decomposition treated in [62], illustrated for three qubits. a) In the first stage 
of the decomposition the operator is broken into two operations controlled on the low qubit, separated by an operation on the 
low qubit controlled on the other qubits. b) The two outer controlled operations are themselves decomposed, resulting in four 
operations on two (in general, n — 1) qubits. In the case of n qubits, these n — 1 qubit unitaries are factored again, until a 
circuit using only one and two qubit gates is obtained. This decomposition is the quantum Shannon decomposition [59], and 
using Cartan involutions one can implement the decomposition given an arbitrary unitary matrix and produce the factors [62] . 



Secondly, we have: 

[fci,mi] € m for k\ in k and ui\ in m. (10) 

finally wc find: 

[nil, "12] G k for toi in and TO2 in m. (11) 

In particular, we note that any subalgebra of m (m itself is not a subalgebra) must be abelian, as all commutators 
must be zero, or they would lie in k, violating the subalgebra property. Hence there may be abelian subalgebras of m 
and we may find such a subalgebra of largest dimension, call it a. It is a fact, that we shall not prove here, that any 
element of m may be written: 

K[hKi = mi, (12) 

where K\ is an element of the subgroup K and h is an element of a maximal abelian subalgebra of to. Furthermore, 
every element of G can be written: 

G\ = K%Mi = exp(&2) expTOi, &2 € k m\ 6 m (13) 

and we may write: 

G\ = exp(fc2)expmi = exp(&2) exp K\hK\, k2,k\ € k m\ em, h G a (14) 
using an elementary property of the exponential map we obtain: 

G\ = exp(fc 2 ) exp K \hK\ = exp{k2)K\exphKi, k 2l ki € k mi € m, h € a (15) 
which we may write: 

G x = K3AK1 (16) 

where K3 — exp(k2)K\ and A = exp(h). To put this in somewhat more physical language, given a Cartan decom- 
position of a unitary group one can find a maximal commuting subalgebra h (i.e. a set of observables that may 
simultaneously take definite values) and a subalgebra K . Then the whole group can be written as a product of 
elements of K and of A, the abelian group corresponding to the maximal abelian subalgebra. For detailed proofs of 
the above assertions we refer the reader to the mathematics literature [63, 64]. 

So far we have only sketched the existence of this decomposition. One further feature of these decompositions 
enables them to be obtained directly for a given element of G and specific K and a. This is the existence of a Cartan 
involution - a self inverse map that fixes the subalgebra k and negates elements of the complement m: 



e(k) = k 

0(m) = — m, 
6(6(x)) = x 



(17) 
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The involution is commutator preserving: 

9([x,y]) = [0{x),0(y)] (18) 

and these maps also extend to the group, with the following useful properties: 

e(K) = K, 0(M) = M -1 , Q(KM) = KM- 1 (19) 

Perhaps the most straightforward example is the involution for the orthogonal subalgebra, where: 

0(x) = -x T , Q(X) = X* (20) 

and the corresponding involution for the case where K is the subalgebra of block diagonal matrices - for example 
where p = q = 2™~ 1 : 

0(x)=( I T J )4 Ip a P J ) (21) 



Irxr J \ Irxr r 

where I pxp is the p x p identity matrix. Similarly 

@ w = (V -L) x (V -L) ^ 

Given these involutions and their properties 19 it is straightforward to obtain M 2 : 

M 2 = (MK-^KM = <d(M^K^)KM = 9(G t )G (23) 

From this point all other factors may be obtained, and recursive application of a sequence of alternating decompositions 
can obtain the quantum shannon decomposition [60, 62]. 

The utility of this decomposition is many-fold. Firstly, for two qubit circuits the accidental isomorphism SU (2) x 
SU(2) — 5*0(4) means that the Cartan decomposition over the orthogonal subalgebra can be used to obtain circuits 
in which the K factors are simply one qubit rotations [65, 66]. Secondly, it was recently shown that the quantum 
Shannon decomposition could be realized by using a pair of Cartan decompositions of the Lie algebra of the Unitary 
group [62]. This insight provided a simple constructive method for obtaining circuits given a unitary operator. This 
method has also been implemented in the python programming language and extensively tested on random unitary 
matrices on up to 8 qubits. The method is suprisingly fast - a fit to timing data gives a rule of thumb timing of 
T = 2 2 - 5 ( g-4 - 5 ) seconds for a q qubit unitary, corresponding to about a minute for a seven qubit unitary operator (a 
128 x 128 matrix). 



III. QUANTUM CHEMISTRY - THE CI METHOD 

The starting point for problems in full-CI calculations is the electronic Hamiltonian expressed in a basis. We briefly 
review the scalable second quantized approach, referred to as the Direct mapping in [25] , which is covered extensively 
in the chapter in this volume by Aspuru-Guzik, before moving on to describe the Compact mapping of [25]. The 
Compact mapping is based on a treatment of the CI matrix, and so we describe the construction of that matrix and 
the method for simulation arising from the use of the CI matrix, we close the section by comparing the direct and 
compact mappings for small numbers of qubits. 



A. Second quantization - the Direct Mapping 

The Hamiltonian for the full quantum treatment of a set of nuclei and electrons is: 

jjmol = f e + f z + Vzz(L pq ) + V ee ( nj ) + V eZ {R pi ) (24) 

where T e and Tz are the electron and nuclear kinetic energies, respectively, V e z(R P i) and V ee (jij) are the electron- 
nuclear and electron-electron coulomb interactions, and Vzz(L pq ) is the nuclear-nuclear coulomb interaction. If we 
invoke the Born-Oppcnheimer approximation and neglect Tz and treat Vzz classically then we obtain: 



N ry N 

If : ■ --- ;Vv! 

= 1 i,L >JJ i>j 



1 lY.^-Y.^ + Y.h (^) 
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which is the Hamiltonian for the electronic structure problem. 

To attack this problem via quantum computation, one requires a mapping of the time evolution operator into a set 
of unitary operations performed on one or two qubits (a quantum circuit). In the direct mapping this is accomplished 
by mapping the states of qubits into the occupation number basis. One therefore requires the same number of qubits 
as spin-orbitals, and a qubit in state one indicates that the corresponding spin-orbital is occupied. The Hamiltonian 
that acts on the fcrmionic states is: 



H 



p,q 



- E 

p,q,r,s 



hpqrs&pQ j qQ>rQ j s •> 



(26) 



where a p is the creation operator for spin-orbital p and a p is the annihilation operator for spin-orbital p. 
The operator a,j acts on basis vectors as follows: 



dj \lo,...,lj-l,l,lj+l,...,lq-l) — (— 1) 

Oj \lpo, . . . , lj-i, 0, ■ . ■ ,l q -i) = 0; 



|Zq, ■ ■ ■ , lj-1,0, lj + 1, • ■ ■ , lq-l) , 



where is the Hermitian conjugate. In order to realize the Hamiltonian 26 as an operator on the qubits we must 
implement the action of the creation and annihilation operators on the qubit states. We do so using the Jordan 
Wigner transformation following [67, 68] . The creation operator is written in terms of Pauli matrices as a} qubit = \ (a x 
- ia y ) = o~ , we can completely express the fermionic annihilation and creation operators in terms of Pauli matrices: 



(27) 
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(28) 



This correspondence is the Jordan- Wigner transformation. The use of this method imposes a number of extra gate 
applications required that scales as 0(q). As the general Hamiltonian 26 has order q terms the use of this method 
imposes a number of gates scaling at least as q 5 , where q is the number of qubits (equivalently of spin-orbitals) . Because 
the terms in the Hamiltonian 26 do not commute one may not simply separate the terms into gates. However, one 
may make use of the Trotter-Suzuki decompositions for small timesteps in order to obtain a circuit with polynomial 
number of gates. The number of gates in such decompositions scales as the number of terms in the Hamiltonian, with 
a prefactor depending on the details of the method, which in turn set the accuracy with which the circuit reproduces 
the time evolution operator. The generic scaling with number of qubits is therefore q 5 for these scalable methods. 



B. FCI - the Compact mapping 

The alternative to second quantized is, of course, the first quantized approach in which all symmetry properties are 
retained in the wavefunction. We take as our starting point a set of 2K spin-orbitals (single electron basis functions), 
for a system of N electrons. From these basis functions we may form 



2K 
N 



(29) 



antisymmetrized TV electron basis functions, or Slater determinants. It is convenient to classify these with reference to 
a Hartree-Fock basis state - the Slater determinant |^o) formed from the N lowest energy spin-orbitals. One may label 
states by their differences from the Hartree-Fock state. The single excitations are labelled \^a), where spin-orbital a 
in |^o) is replaced by spin-orbital r. Similarly, doubles and triples are labelled l^alc) an d so on - The number 

of n-tuple excitations is: 



N\ (2K - AT 
n ) V n 



(30) 



summing over n gives the total number of Slater determinants: 

m t(1)( 2K : N ) <-> 



N ) ' \ n J \ n 

n=0 
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a result obtained from the Vandcrmondc convolution and the symmetry property of binomial coefficients. 

Given this basis of determinants one wishes to obtain the Hamiltonian matrix elements within this basis - the CI 
matrix. To obtain the matrix elements of the Hamiltonian one uses Slater's rules [69-71]. Given a matrix element 
(A\H\B), where A and B are configurations (lists of spin-orbital labels), we may compute the CI matrix elements as 
follows: 

1. Diagonal elements. If A = B so that the determinants are identical 

n n / \ 

(A\ H\A)=J2 (P\ h \ P ) + £(<P9lM> - <pg|«p> J (32) 

P P<Q 



2. If the configurations A and B differ in one spin-orbital, p in A and q in B 

(A\ H \B) = ( P \ h\q)+J2 (<PJM) - (Pl\k) ) (33) 



3. If the configurations A and B differ in two spin-orbitals, pq in A and rs in B 

(A\H\B) = (pq\rs)-(pq\sr) (34) 

4. If the configurations A and B differ in more than two spin-orbitals 

(A\ H\B)=0 (35) 

where we have defined the one-electron integrals: 

{plhlq) =~^ / WV^Wdx (36) 

and the two electron integrals: 

(pq\rs) = I 0*(xi)0*(x 2 ), r0 r (xi)(/) s (x2)dxidx2. (37) 

J l x l — x 2| 

At this stage quantum simulation of FCI via the compact mapping comes into focus. The compact mapping is a 
mapping of the logical states of the qubits onto labels of Slater determinants. The matrix eigenvalue problem we wish 
to solve is defined by the CI matrix, and we implement time evolution according to this Hamiltonian by, for example, 
using a Cartan decomposition of the relevant time evolution operator. Ironically, while FCI is typically the classical 
method that uses the largest number of states for a given problem, in the quantum case the compact mapping uses 
a smaller Hilbert space dimension than the direct mapping of the second quantized approach. The reason for this is 
that the second quantized approach maps the entire Fock space of the system into the Hilbert space of the qubits, 
whereas the compact mapping only maps the TV electron sector into the quantum computer. 

The procedure outlined above for the compact mapping is not scalable. Although the input data to the CI matrix is 
the same set of one and two electron integrals used to define the second quantized Hamiltonian, the compact mapping 
does not respect the tensor product structure of the qubit Hilbert space, and hence there is no known procedure 
analogous to the Trotter-Suzuki decomposition that enables the efficient construction of circuits with a number of 
gates rising polynomially in the number of spin-orbitals. In addition, while for the second quantized formalism the 
one and two electron integrals appear directly in the circuit as rotation angles for quantum gates acting on a small 
number of qubits, at present there is no known way of constructing the time evolution operator corresponding to the 
CI matrix on the quantum computer without computing it classically first. 

The compact mapping remains of interest for a number of reasons, however. If we consider problems in which we 
have 2K spin-orbitals, and therefore the second quantized approach requires 2K qubits (for up to 2K electrons), the 
compact mapping requires: 



c(k,N) = 



log 2 



(38) 



which is less than 2K for all N < 2K, being the nearest whole number of qubits into which the N particle sector 
of the full fock space of up to 2K particles will fit. Hence the compact mapping will always enable the mapping of 
larger problems than the direct mapping into the qubit Hilbert space. 



9 



This advantage is only of interest if one can in fact take advantage of the smaller qubit requirements by implementing 
the time evolution operator. Evidently, the non-scalable nature of the compact mapping makes it asymptotically 
uninteresting. However, the issue for small examples is not quite so clear cut. Suppose the number of gates required 
to implement the evolution operator which yields the U is A(2K) 5 , where A is a constant. In order to obtain p bits 
of precision in the final answer one must run at least one phase estimation procedure with 2 P A(2K) 5 gates. The 
corresponding numbers for performing the Cartan involution are much worse for the first stage of iterative phase 
estimation, where a Cartan involution requires 



22c(k.N) _ 



(39) 



gates, as well as a comparable amount of classical computation to obtain the circuit. However, each stage of the 
iterative phase estimation requires only the same number of gates, as one simply computes new parameters in the 
Cartan decomposition for each evolution operator U 2 . Hence the total gate count comparison is roughly: 



Gc 



ompact 



(40) 



compared with 



Gdi: 



2 P A(2K) 5 



(41) 



The crossover point at which the advantage in number of gates of the direct mapping outweighs the advantage in 
qubits of the compact mapping represents a significant milestone for quantum simulation - the point at which the 
advantage of scalable methods is realized in practice. We plot G compac t and GDirect in Figure 3 for A = 1, for 5 and 10 
electrons and qubit numbers between 5 and 20. These estimates indicate that the crossover point between Compact 
and Direct mappings is between 10 and 20 qubits. 
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FIG. 3: Left: A comparison of gate counts, as estimated by equations (40) and (41) (with A = 1) for 5 electrons and numbers 
of qubits between 5 and 10. Solid line: GDirect- Dashed line, Gcompact- This indicates that for fewer than 5 qubits the Compact 
mapping is preferred over the scalable Direct mapping. Right: A comparison of gate counts, as estimated by equations 40 
and 41 (with A = 1) for numbers of qubits between 5 and 10. Solid line: GDirect, 5 electrons. Dashed line, Gcompact, 5 electrons. 
Dotted line, Gcompact, 10 electrons. This indicates that the crossover between the Compact mapping and the scalable Direct 
mapping lies between 10 and 20 qubits. 



IV. A SELECTION OF HISTORICAL CALCULATIONS IN QUANTUM CHEMISTRY 

Here we present a sequence of calculations spanning from the 1930's to the 1960's. Our intent here is not to give a 
historically faithful account of the development of quantum chemistry. We merely choose a few mileposts and evaluate 
the cost in qubits of the Compact mapping defined above of those calculations on quantum hardware. The choice of 
these calculations was strongly influenced by the presentation of [72] . For translations of earlier foundational papers 
in quantum chemistry we refer the reader to Hinne Hettema's volume [73], which also includes a fascinating preface 
discussing the historical context of early work on quantum chemistry in both Germany and the United States. 
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A. 



The thirties and forties 



Perhaps the first exact ab initio calculation was performed by Coulson on a Brunsviga mechanical calculator for the 
molecular ion H$ [106] . Early calculations of the ground state energy of the Hydrogen molecule were conducted in the 
1930's using the variational principle [74]. The first studies by Heitler and London were superceded by the methods 
of Hyelleras in which a variational Ansatz is constructed from symmetry considerations alone. The calculation of the 
variational energy in these calculations naturally involved integrals over the Ansatz functions. 



In 1950 Boys introduced the use of Gaussian functions [75]. In the same year Boys gave a presentation of the 
CI method and presented a CI calculation of the ground state energy of Beryllium [76] . Interestingly, Boys first CI 
calculation used Slater-type LCAO orbitals rather than Gaussians because of the relatively simpler integrals which 
occur in this atomic problem. He used three s and one p single-electron functions and performed the calculation using 
six and ten determinants formed from these functions. Here we see perhaps the first application of CI, corresponding 
to a three or four qubit problem. 

In 1952 Taylor and Parr applied the CI method to the Helium atom, a system which they note is "no longer an 
interesting problem per se" because of the accuracy of Hycllcraas' variational solution obtained in 1930. However 
they note that in spite of the melancholy prospects for progress of the CI method Helium is the first non-trivial 
example. Their calculation involved four configurations lsls', 2p 2 , 3d 2 , 4/ 2 and hence corresponds to a two qubit 
calculation. Their paper reproduces 88% of the correlation energy of Helium [77]. Three years later, Shull and Lowdin 
reconsidered the Helium atom using natural orbitals [78]. They obtained the CI solution using six configurations from 
three orbitals (Is, 2s, 3s). The results were then analyzed in terms of natural orbitals, and they showed that similar 
results could be obtained with only three configurations. By doing so they reduced a three-qubit (six configuration) 
calculation to a two-qubit (three configuration) calculation. 

In 1956 CI calculations were reported on BH with 7 orbitals, 23 co-detors, and 3 internuclear distances, H 2 with 
8 orbitals, 96 co-detors, at seven nuclear configurations. These results were obtained by Boys using the EDSAC 
computer at Cambridge [41]. The BH calculations are five-qubit problems, and the H 2 calculations are seven 
qubit problems. The following year Miller, Friedman, Hurst and Matsen continue the extension of CI to molecular 
problems with calculations of LiH and BeH + [79]. Their calculations used 20 terms corresponding to all possible 
singlet states of Is, 2s, 2p a Slater orbitals on the metal atoms and Is on the Hydrogen. They performed six, ten and 
twenty determinant calculations on IBM CPC and IBM 650 machines, these are three, four and five qubit problems 
respectively. In the same year, Pekeris gave a new method for the determination of the ground state energy of 
two-electron atoms [80]. 

In a largely philosophical paper published in 1959 Mulliken and Roothaan reflected on the status of exact quantum 
chemistry noting that "it so happens that the more difficult problems of chemistry and molecular physics are of a 
high order of complexity in terms of the mathematical analysis and the computational efforts required for quantitative 
results. However, and this is the point of our paper, those who are working in this field are now catching up, and we 
believe that the results of theoretical computations are going to compete more and more strongly with experiment 
from now on" [81]. 



In 1960 Watson published a CI study of Be involving 37 configurations constructed from a basis of STO's- a six qubit 
problem [82]. Boys and Cook also lay out the progress made in developing the CI method over the previous decade 
in [95]. The following year Sinanoglu published perhaps the first in depth study of electron correlation, introducing 
the idea that the correlation energy would be well approximated by the sum of pair correlation energies [96] . 

In 1963 Barnctt at MIT reports on the development of a FORTRAN code for the IBM 709/90 which implements 
SCF and CI calculations on basis sets of size < 50 although no actual example calculations apart from integral 
evaluations are reported in this work [97]. In the same year Hagstrom and Shull report on the efficacy of natural 
orbitals for H 2 and H 2 . They report that calculations using 21 and 33 configurations may be replaced by natural 
orbital bases of size 8 and 15 respectively [84]. For our purposes this is of particular interest as it converts five and 
six qubit calculations into three and four qubit calculations. Natural orbital calculations for HeH and Li 2 followed in 
1966 providing seven configuration (or three qubit) calculations [85]. These calculations are of interest as the natural 
orbital wavefunctions provide a class of ansatz states whose complexity of preparation can be investigated. 



B. 



The fifties 



C. 



The sixties 
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Year 


Calculation 


Citation 


Number of qubits 


1933 


H 2 


[74] 


1 


1950 


Be 


[76] 


3,4 


1952 


He 


[77] 


2 


1955 


He 


[78] 


2,3 


1956 


BH 


[41] 


5 


1956 


H 2 


[41] 


7 


1957 


LiH 


[79] 


3,4,5 


1957 


BeH+ 


[79] 


3,4,5 


1960 


Be 


[82] 


6 


1960 


CH 2 


[83] 


19 


1963 


H 2 


[84] 


3,4,5,6 


1966 


HeH 


[85] 


3 


1966 


Li 2 


[85] 


3 


1967 


H 2 


[86] 


10 


1967 


H 2 


[87] 


24 


1967 


H 2 


[88, 89] 


38,39 


1968 


H 2 


[90] 


39,46 


1968 


Be 


[91] 


11 


1969 


Li, Be+, B ++ 


[92] 


9,10 


1969 


BH, FH 


[93] 


12, 14 


1970 


H 2 


[94] 


23 



TABLE I: A sequence of quantum chemical problems and their minimal qubit requirements. These calculations bridge the 
period of development of CI from the first variational calculations of the hydrogen molecule energy to the advent of the 
GAUSSIAN- 70 program. The qubit costs are associated with the Compact mapping of full-CI calculations on the system 
studied in the cited publication, although the methods used originally may not be full-CI. 

In the same year (1968) Bunge published a new calculation of the ground state energy of the Be atom, using a 180 
term CI with 1492 Slater determinants - an eleven qubit calculation [91]. A study in 1969 of the irir* states of the 
ethylene molecule compared HF and CI calculations [98]. An improved CI algorithm was applied to three-electron 
problems (Li, Be + and B ++ ) using 7 to 9 Slater orbitals, (9 and 10 qubit problems respectively) [92]. Calculations on 
i?3 were also performed with 9 Slater orbitals, another 10 qubit problem. The publication by Bender and Davidson 
in 1969 of a study of the first row diatomic hydrides using very large CI provides a rich set of examples [93]. Two 
example calculations are of BH, a twelve qubit calculation, and of FH, a fourteen qubit calculation. 

In 1970, a new study of water with a 14 function basis set, corresponding to a 23 qubit calculation, comparing 
STO and contracted Gaussian basis sets [94]. In the same year one of the first calculations using the GAUSSIAN- 70 
code appeared - treating 30 species with STO-3G basis sets in which Slater type functions are approximated by three 
Gaussians [99]. With this paper and the advent of codes such as GAUSSIAN which achieved wide usage in the 
theoretical quantum chemistry community full CI calculations became common. The Pople classification of basis sets 
makes the determination of a number of qubits for any large CI calculation straightforward, and the number of qubits 
is typically 20 — 50 for calculations which historically lie between 1970 and the present. These basis sets and relevant 
citations are readily available from, e. g. [100]. 

In 1960 several papers applying variational and LCAO-SCF methods to molecules appeared. Mclean reports LCAO- 
SCF on C2H2 and CO2 with molecular orbitals constructed from (for C2H2) Is STO's on each H and Is, 2s, 2pa 
and 2p7r STO's on each C. For C0 2 MO's are formed from Is, 2s, 2pa and 2pn STO's on each O and Is, 2s, 2pa 
and 2pn on the C [102]. In the same year there was also a calculation of the excited state of H 2 by a variational 
method [103]. 

In 1960 Nesbet reports an LCAO Hartree-Fock calculation on HF using 9er and 5tt orbitals using the IBM 704 
computer [104]. Foster and Boys report a study of the three lowest lying states of 6 different bond angles for the CH 2 
radical using a basis of six exponential functions on the carbon and one on each hydrogen [83]. A CI calculation in 
this basis would represent a 19 qubit problem. 

Much work in the previous century was devoted to approaching the Hartree-Fock limit using ever better SCF 
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Calculation 


Method 


Size of basis 


Energy (a.u.) 


Number of qubits 


Reeves and Boys [101] 


MC STO 




-75.776 




McWeeny and Ohno [86] 


MC STO 


7 


-75.761 


10 


Moccia 


OC STO 


28 


-75.992 


36 


Harrison [87] 


MC GLF 


14 


-76.002 


24 


Moskowitz and Harrison [89] 


MC GTO 


36 


-76.034 


39 


Ritchie and King [88] 


MC CGF 


38 


-76.034 


38 


[90] 


MC GLF 


32 


-76.044 


39 


[90] 


MC GTO 


56 


-76.002 


46 



TABLE II: SCF calculations of the water molecule ground state energy in 1967 and 1968. These calculations, aiming to approach 
the Hartree-Fock limit for water define large molecular orbital bases. The qubit numbers refer to the minimal requirements 
to perform a full-CI in the molecular orbital basis defined by the published SCF calculation. Unsuprisingly, SCF calculations 
which approach the Hartree-Fock limit define full-CI problems requiring several tens of qubits. 



methods. For the water molecule, an example of this work is given by [90]. In this work, progress to date (prior 
to 1968) for the SCF treatment of water is summarized in Table IX, which we reproduce in Table II. These SCF 
calculations provide an interesting reference: repeating these calculations using CI yields more challenging, but still 
tractable problems. The data in Table II give us qubit requirements in the range 10 — 46 



V. BOYS' 1950 CALCULATION FOR BE 



In this paper Boys constructs a 12 x 12 CF Hamiltonian over a set of configuration state functions (CSF's) which 
are referred to as "co-detors" - spin-adapted linear combinations of determinants. There are four sets of functions 
sA, sB, sC, and pA, where the s and p notation indicates that the functions in these sets arc eigenvalues of the 
square of total angular momentum with quantum number I = (s) and I = 1 (p). The letters A, B, C denote 
different radial dependence, and for each radial dependence and total orbital angular momentum there is one function 
for each possible combination of the eigenvalues of L z and S z . Because these are atomic orbitals (single electron 
eigenfunctions) the total spin eigenvalue is 1/2. Writing out the elements of each set we obtain: 

sA = {sA(0, -1/2), sA(0, 1/2)} 
S B = {sB(0,-l/2),sB(0,l/2)} 

sC = {sC(0, -1/2), sC(0, 1/2)} [ ' 

P A = {pA(-l, -l/2),pA(-l, 1/2), pA(0, -1/2),M(0, l/2),pA(l, -l/2),pA(l, 1/2)} 

Hence there are 12 atomic orbitals here. However, Boys uses linear combinations of functions in this basis which are 
simultaneous zero eigenvalues of I? and S 2 . He does so by taking the four electrons and placing them in two pairs, 
each pair being in a spin and orbital singlet. 

The orbital coupling can use the special case of the Clcbsch-Gordan coefficients relevant to the singlet state of total 
orbital angular momentum, (Zi, l 2 , mi, m 2 \J, M), where li,l 2 ,mi, M 2 are the angular momentum quantum numbers 
labeling the orbitals and J and M are the total angular moment quantum numbers. By the recoupling rules for 
angular momentum such a state can only be produced from a pair of angular momentum states with equal values of 
I: 

(l,l,m 1 ,m 2 \0,0)=S mu _ m2 ^=,. (43) 

This implies that the combination of orbital products is a uniform superposition of pairs with opposite z components 
with alternating signs. 



1 1 (_\\l-mi 

l°>°) = H H 5 mi ^ m2 —j==-\l,l,mi,m 2 ) 



mi— — I m2 — — l 



= £ \l,l,mi,-mi) 

mi=—l 



(44) 
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These are denoted: 



(sAsAYSisBsBYS, 
(sA8A) 1 S(sBsC) 1 S, 
(sAsA^SipApA^S, 
(sAsB) 1 S(pApA) 1 S 



(45) 



These are then antisymmetrized by applying the operator: 



S = signTr^L 



4! 



(46) 



Resulting in the normalized wavefunctions given in Eqn 4 in [76]. 

The elementary integrals (overlaps and expectation values for the one-electron orbitals) whose radial parts are given 
by: 



(r\A,l,a) = r A+l e- ar 



These radial functions are not orthogonal. Their overlaps are given by: 

The values of these parameters for the functions sA, sB, sC and pA are given in table III. 



(47) 



(48) 



Orbital 


\A,l,a) 


(r\A,l,a) 


A 


I 


a 


{0A\l,rn) 


m 


sA 


10,0,4} 


e~ 4r 








4 







sB 


|1,0,1) 


re~ r 


1 





1 


rr_ 





sC 


10,0,3} 


e~ 3r 








3 


V 4-7T 





pA 


|0,1,1) 


re~ r 





1 


1 




1 


pA 


10,1,1} 


re~ r 





1 


1 


\ / 7T^- cos 9 

y 





pA 


10,1,1} 


re~ r 





1 


1 


y c" 71 " 


-1 



TABLE III: Hydrogenic orbitals for the calculation of [76] 



Orbital 


\A,l,a) 


10,0,4} 


11,0,1} 


|0,0,3) 


sA 


(0,0, 4| 


(1/16) 2 


3!/5 4 


2!/7 3 


sB 


(1,0, 1| 




4!/2 5 


3!/4 4 


sC 


(0,0, 3| 






2!/6 3 



TABLE IV: Overlap integrals for the calculation of [76] 



The three p orbitals are orthogonal to the s orbitals and to each other by virtue of the orthogonality of the angular 
parts of the wavefunction. The three s orbitals are not orthogonal, but an orthogonal basis may be constructed from 
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them by the Gram-Schmidt procedure. The corresponding orthonormal orbitals are: 

1 



\<h) 



V(0,0,4|0,0,4) 

(0,0,4|1,0,1 



2> = |1,0,1) 
1 



|0,0,4) = 16|0,0,4) 
|0, 0, 4) 



\h) = 



(0,0,4|0,0,4) 

3 5 4 7 3 - 2 17 



1 



yp^T L 54 255473 - 32l773 /5 4 ^jmm 73 

1 5 4 7 3 - 2 17 

^/(^) 255473-3x2"73/5- 

4 |o,i,i) 



|0,0,4) 



(49) 



|1,0,1) + |0,0,3) 



Where the coefficients in eqn (18) in [76] can now be seen to arise as follows: 

1 



(0,0,4|1,0,1) 



3! 



x/W'alV's) 



V ^V^ 2 )(0,0,4|0,0,4) 

3 5 4 7 3 - 2 17 1 1 

54 255 4 73 - 32i 7 73/5 4 ~ 73 

1 S 4 ? 3 - 2i7 

7P^y255473 - 3 x 2i773/54 

1 



After all that, the six by six CI matrix actually computed by Boys was: 



1.173302451 

2.883508103 

-69.59716965 

0.5968868832 
47.60738095 

1.154700538 



(50) 



/ 0. 0.3967 0.2488 -0.9416 -0.1361 0. \ 

0.3967 17.0628 -0.1257 -0.5601 0. 0. 

0.2488 -0.1257 5.0387 0.1499 -0.0099 0.0196 

-0.9416 -0.5601 0.1499 8.7275 0. -0.1361 

-0.1361 0. -0.0099 0. 0.3677 -0.9502 

V 0. 0. 0.0196 -0.1361 -0.9502 9.1242 / 



(51) 



Given the CI matrix and established a compact mapping for the co-detors we can consider a non-scalable imple- 
mentation of Boys original calculation. Our starting point is the 6x6 version of the CI matrix. As noted above, the 
task is to decompose powers of the evolution operator corresponding to this Hcrmitian matrix into a quantum circuit. 
This may be accomplished by a Cartan decomposition that treats the evolution operator as an action on a qubit 
and a qutrit. While it is conventional to decompose unitaries over qubits, controlled quantum systems employing 
larger numbers of states than two have also been utilized in a number of experimental realizations [105]. A more 
conventional approach would utilize three qubits in the state register with two three qubit states uninvolved in the 
calculation. The optimal three-qubit circuit requires 44 CNOT gates and 63 one qubit gates [62]. 



VI. CONCLUSION 



Elucidation of small cases of quantum chemical calculations from the past literature provides an interesting exper- 
imental road map leading from calculations which quantum computers can perform now, to calculations beyond the 
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Boys Designation 


Qubit state 


A(sAsAyS(sBsBys 


|0> 


AisBsB^SisCsCyS 


|1> 


AisAsA^SisBsCyS 


|2> 


AisAsCySisBsB^S 


|3> 


A(sAsA) 1 S(pApA) 1 5 


|4) 


AisAsCySipApA^S 


|5> 


A(sAsAys(sCsCys 


|6) 


A(sAsByS(sCsCyS 


|7> 


A(sBsByS(pApAyS 


|8> 


A(sAsB) 1 SipApA) 1 S 


|9> 



TABLE V: Compact mapping of Boys basis for Beryllium 





|0> 


|1> 


|2> 


|3> 


|4> 


|5> 


|6> 


|7> 


|8> 


|9> 




-14.4577 


0.3967 


0.2488 


-0.9416 


-0.1361 


0. 


0.0541 


0.0986 


-0.1361 


0.0276 


11) 


0.3967 


2.6051 


-0.1257 


-0.5601 


0. 


0. 


0.0220 


-0.1958 


-0.0317 





|2> 


0.2488 


-0.1257 


-9.419 


0.1499 


-0.0099 


0.0196 


0.3534 


1.0192 





0.0070 


|3) 


-0.9416 


-0.5601 


0.1499 


-5.7302 


0. 


-0.1361 


0.0269 


-0.8085 





-0.1780 


|4> 


-0.1361 


0. 


-0.0099 


0. 


-14.09 


-0.9502 


-0.0132 





0.0220 


-0.2816 


15) 


0. 


0. 


0.0196 


-0.1361 


-0.9502 


-5.3335 


0.0215 


0.0070 


-0.0269 


0.0407 


|6> 


0.0541 


0.0220 


0.3534 


0.0269 


-0.0132 


0.0215 


-3.7904 


1.3259 








|7) 


0.0986 


-0.1958 


1.0192 


-0.8085 





0.0070 


1.3259 


1.5002 





-0.0946 


18) 


-0.1361 


-0.0317 








0.0220 


-0.0269 








-4.1937 


0.0032 


|9> 


0.0276 





0.0070 


-0.1780 


-0.2816 


0.0407 





0.0946 


0.0032 


-10.1121 



TABLE VI: CI matrix for Boys basis for Beryllium 



capabilities of comparable classical methods. Scalable methods for the calculation of ground state properties using 
quantum computers have been studied extensively in the last few years. The known scalable techniques for implement- 
ing time evolution operators of fcrmionic systems, comprising the Jordan Wigncr transformation and Trotter-Suzuki 
decompositions, provide polynomial performance. However, they are most likely not yet optimized and further study 
of these scalable approaches is likely to be fruitful. At the other extreme, the approach detailed here considers the 
smallest calculations where non-scalable approaches that minimize the required number of qubits are used. These 
allow for early experimental implementation, the purpose of which is to demonstrate sufficient quantum control and 
coherence to obtain results of chemical accuracy for small systems. As we have argued here, there is a crossover 
point between scalable and non-scalable implementations which will likely result in two experimental frontiers. At 
one frontier as much chemistry as possible will be crammed into the available qubit Hilbert space as possible, whether 
the approach used is scalable or not. At the other frontier we will be proceeding with the development of the scalable 
methods that we know will carry us beyond model systems and ultimately beyond the capabilities of any conceivable 
classical computing device. 
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